function varargout = PlotNB_Gaps(ax,sch,cfg)
% PlotNB_Gaps(ax,sch)
% PlotNB_Gaps(ax,sch,cfg)
% HR = PlotNB_Gaps(ax,sch,...)
% [HR,FA] = PlotNB_Gaps(ax,sch,...)
% [HR,FA,GD] = PlotNB_Gaps(ax,sch,...)
% [HR,FA,GD,GDI] = PlotNB_Gaps(ax,sch,...)
% [HR,FA,GD,GDI,S] = PlotNB_Gaps(ax,sch,...)
% [HR,FA,GD,GDI,S,cHR] = PlotNB_Gaps(ax,sch,...)
% 
% ax = handle to axes object
% sch = schedule structure
% 
% cfg fields:
%   .legend = true or false
% 
% 
% DJS (c) 2011

if ~exist('cfg','var') || isempty(cfg)
    cfg = struct('legend',true);
end


if isempty(ax), ax = gca; end

h = 2;      % hits
% m = 1;      % misses
f = 3;      % false alarm

rv = sch.response_vals;

% find trial types
t = strfind(sch.readparams,'trial_type');
for i = 1:length(t)
    if ~isempty(t{i}), xtt = i+1; break; end
end
catch_trials = rv(:,xtt) == 1;

% find response codes
t = strfind(sch.readparams,'response_code');
for i = 1:length(t)
    if ~isempty(t{i}), xrc = i+1; break; end
end

% find HP frequencies of noise stimuli
t = strfind(sch.readparams,'NoiseHP');
for i = 1:length(t)
    if ~isempty(t{i}), xhp = i+1; break; end
end
uhp = unique(rv(~catch_trials,xhp)); uhp(uhp==0) = [];

% find Gap Durations
t = strfind(sch.readparams,'stim_duration');
for i = 1:length(t)
    if ~isempty(t{i}), xsd = i+1; break;    end
end
usd = unique(rv(~catch_trials,xsd));

% compute hit rate for each Frequency & DB level
hr = zeros(length(usd),length(uhp));
for i = 1:length(uhp)
    for j = 1:length(usd)
        ind   = rv(:,xhp) == uhp(i) & rv(:,xsd) == usd(j);
        hind  = ismember(rv(ind,xrc),h);
        hr(j,i) = sum(hind)/sum(ind);
    end
end
hr(isnan(hr)) = 0;
% compute false-alarm rate
fa = sum(rv(catch_trials,xrc) == f) / sum(catch_trials);

% plot data
cla(ax);

% plot(ax,usd,hr,'LineWidth',2);
plot(ax,usd,hr,'o-','MarkerSize',3);

hold(ax,'on');

% compute corrected hit rate and estimate sensitivity
S = nan(size(hr,2)); cHR = nan(sum(usd>0),size(hr,2));
if size(hr,1) > 2 && any(any(hr))
    for i = 1:size(hr,2)
        cHR(:,i) = hr(usd>0,i) - hr(usd>0,i) * fa; % corrected hit rate (Heffner & Heffner, 1988)
        S(i) = interp1q(cHR(:,i),usd(usd>0),0.5); % estimate sensitivity
    end
end

ylim(ax,[0 1]);

set(ax,'xscale','log');

plot(ax,xlim(ax),[fa fa],':r');

hold(ax,'off');

xlabel(ax,'Gap Duration (ms)');
ylabel(ax,'Hit Rate');

if cfg.legend
    h = legend(ax,num2str(uhp/1000,'%1.0fkHz'),'Location','SouthEast');
    set(h,'Box','off');
    set(h,'Color','none');
end

varargout{1} = hr;
varargout{2} = fa;
varargout{3} = usd;
varargout{4} = xsd;
varargout{5} = S;
varargout{6} = cHR;





